#Appendix: Quick. But Impactful? United Nations Quick Impact Projects and Violence against Civilians in Civil War
rm(list=ls())
setwd("YOUR WORKING DIRECTORY HERE")

library(tidyverse)
library(stargazer)
library(readxl)
library(MASS)
library(ggeffects)
library(jtools)
library(marginaleffects)
library(ggpubr)
library(MatchIt)

#open data 
merge <- read.csv("qip_civilians_replication_tables_figures.csv", row.names = NULL)

#subset to only observations that had at least one attack on civilians & at least one project 
temp <- merge %>% 
  group_by(admin)%>%
  filter(sd(rebel.violence.against.civilians.lead1, na.rm = TRUE) > 0 & sd(log.cost_qipongoing, na.rm = TRUE) > 0) %>% 
  ungroup()

#TABLE A1: Main models
summary(m1 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.cost_qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m2 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))

stargazer(m1, m2, 
          covariate.labels = c("QIPs-cost (logged)", "QIPs-number (logged)", "Rebel violence against civilians", "Government violence against civilians", "Battles", "Mission duration", "PK troops (logged)", "Number of armed actors"), 
          add.lines = list(c("Fixed Effects", "District", "District")),
          type="html", omit="admin", out="out.doc")

#TABLE A2: Where do QIPs go?
summary(m1 <- glm.nb(qipongoing ~ battles, data=merge))
summary(m2 <- glm.nb(qipongoing ~ log.pktroops, data=merge))
summary(m3 <- glm.nb(qipongoing ~ pkbase, data=merge))

stargazer(m1, m2, m3,  
          covariate.labels = c("Battles", "PK troops (logged)", "PK bases"), 
          add.lines = list(c("Fixed Effects", "None", "None", "None", "None")),
          type="html", out="out.doc")

#TABLE A3: Controlling for the number of peacekeeping bases
summary(m1 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.cost_qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m2 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))

stargazer(m1, m2, 
          covariate.labels = c("QIPs-cost (logged)", "QIPs-number (logged)", "Rebel violence against civilians", "Government violence against civilians", "Battles", "Mission duration", "PK troops (logged)", "Number of armed actors", "Number of PK bases"), 
          add.lines = list(c("Fixed Effects", "District", "District")),
          type="html", omit="admin", out="out.doc")

#TABLE A4: Violence against civilians lags & cost of QIPs
summary(m1 <- glm.nb(cost_qipongoing ~ vac.lag1 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m2 <- glm.nb(cost_qipongoing ~ vac.lag2 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m3 <- glm.nb(cost_qipongoing ~ vac.lag3 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m4 <- glm.nb(cost_qipongoing ~ vac.lag4 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m5 <- glm.nb(cost_qipongoing ~ vac.lag5 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m6 <- glm.nb(cost_qipongoing ~ vac.lag6 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))

stargazer(m1, m2, m3, m4, m5, m6, 
          covariate.labels = c("Violence against civilians (1 mo lag)", "Violence against civilians (2 mo lag)", "Violence against civilians (3 mo lag)", "Violence against civilians (4 mo lag)", "Violence against civilians (5 mo lag)", "Violence against civilians (6 mo lag)", "Government violence against civilians", "Battles", "Mission duration", "PK troops (logged)", "Number of armed actors"), 
          add.lines = list(c("Fixed Effects", "District", "District", "District", "District", "District", "District")),
          type="html", omit="admin", out="out.doc")

#TABLE A5: Violence against civilians lags & number of QIPs
summary(m1 <- glm.nb(qipongoing ~ vac.lag1 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m2 <- glm.nb(qipongoing ~ vac.lag2 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m3 <- glm.nb(qipongoing ~ vac.lag3 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m4 <- glm.nb(qipongoing ~ vac.lag4 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m5 <- glm.nb(qipongoing ~ vac.lag5 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m6 <- glm.nb(qipongoing ~ vac.lag6 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))

stargazer(m1, m2, m3, m4, m5, m6, 
          covariate.labels = c("Violence against civilians (1 mo lag)", "Violence against civilians (2 mo lag)", "Violence against civilians (3 mo lag)", "Violence against civilians (4 mo lag)", "Violence against civilians (5 mo lag)", "Violence against civilians (6 mo lag)", "Government violence against civilians", "Battles", "Mission duration", "PK troops (logged)", "Number of armed actors"), 
          add.lines = list(c("Fixed Effects", "District", "District", "District", "District", "District", "District")),
          type="html", omit="admin", out="out.doc")

#TABLE A6: Matching: at bottom of document

#TABLE A7: Alternative leads (violence against civilians at time t+2 and t+3)
summary(m1 <- glm.nb(rebel.violence.against.civilians.lead2 ~ log.cost_qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m2 <- glm.nb(rebel.violence.against.civilians.lead2 ~ log.qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m3 <- glm.nb(rebel.violence.against.civilians.lead3 ~ log.cost_qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m4 <- glm.nb(rebel.violence.against.civilians.lead3 ~ log.qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))

stargazer(m1, m2, m3, m4,
          add.lines = list(c("Fixed Effects", "District", "District", "District", "District")),
          covariate.labels = c("QIPs-cost (logged)", "QIPs-number (logged)", "Rebel violence against civilians", "Government violence against civilians", "Battles", "Mission duration", "PK troops", "Number of armed actors"), 
          type="html", omit="admin", out="out.doc")

#TABLE A8: Alternative measures of baseline violence
summary(m1 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.qipongoing + vac.average.3mo + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m2 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.cost_qipongoing + vac.average.3mo + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m3 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.qipongoing + vac.lag12 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m4 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.cost_qipongoing + vac.lag12 + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))

stargazer(m1, m2, m3, m4,
          add.lines = list(c("Fixed Effects", "District", "District", "District", "District")),
          covariate.labels = c("QIPs-cost (logged)", "QIPs-number (logged)", "3 month average-rebel violence against civilians", "1 year prior-rebel violence against civilians", "Government violence against civilians", "Battles", "Mission duration", "PK troops", "Number of armed actors"), 
          type="html", omit="admin", out="out.doc")

#TABLE A9: Decaying value of QIPs
summary(m1 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.cost_qipongoing_decay + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m2 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.qipongoing_decay + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops +  n_actors + factor(admin), data=temp))

stargazer(m1, m2,
          covariate.labels = c("QIPs-cost (logged)", "QIPs-number (logged)", "Rebel violence against civilians", "Government violence against civilians", "Battles", "Mission duration", "PK troops", "Number of armed actors"), 
          add.lines = list(c("Fixed Effects", "District", "District")),
          type="html", omit="admin", out="out.doc")

#TABLE A10: The relationship between other aid and violence against civilians
temp.otheraid <- merge %>% 
  group_by(admin)%>%
  filter(sd(rebel.violence.against.civilians.lead1, na.rm = TRUE) > 0 & sd(log.other.aid.projects, na.rm = TRUE) > 0) %>% 
  ungroup()

summary(m1 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.other.aid.projects + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp.otheraid))
stargazer(m1,
          covariate.labels = c("Other types of aid (logged)", "Rebel violence against civilians", "Government violence against civilians", "Battles", "Mission duration", "PK troops", "Number of armed actors"), 
          add.lines = list(c("Fixed Effects", "District")),
          type="html", omit="admin", out="out.doc")

#TABLE A11: Country fixed effects
summary(m1 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin) + factor(country), data=temp))
summary(m2 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.cost_qipongoing + rebel.violence.against.civilians + battles + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin) + factor(country), data=temp))

stargazer(m1, m2,
          covariate.labels = c("QIPs-cost (logged)", "QIPs-number (logged)", "Rebel violence against civilians", "Government violence against civilians", "Battles", "Mission duration", "PK troops", "Number of armed actors"), 
          add.lines = list(c("Fixed Effects", "District & Country", "District & Country")),
          type="html", omit=c("admin", "country"), out="out.doc")


#TABLE A6: Matching
cov <- c('rebel.violence.against.civilians', 'gov.violence.against.civilians', 'battles', 'mission.duration', 'log.pktroops', 'n_actors')
temp %>%
  group_by(qipongoing.binary) %>%
  dplyr::select(one_of(cov)) %>%
  summarise_all(funs(mean(., na.rm = T)))
m_ps <- glm(qipongoing.binary ~ rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors,
            family = binomial(), data = temp)

prs_df <- data.frame(pr_score = predict(m_ps, type = "response"),
                     qipongoing.binary = m_ps$model$qipongoing.binary)

#examining region of common support
labs <- paste("Ongoing QIPs:", c("Yes", "No"))
prs_df %>%
  mutate(qipongoing.binary = ifelse(qipongoing.binary == 1, labs[1], labs[2])) %>%
  ggplot(aes(x = pr_score)) +
  geom_histogram(color = "white") +
  facet_wrap(~qipongoing.binary) +
  xlab("Probability of receiving QIP") +
  theme_bw()

#executing matching algorithm
temp_nomiss <- temp %>%  # MatchIt does not allow missing values
  dplyr::select(rebel.violence.against.civilians.lead1, qipongoing.binary, log.cost_qipongoing, log.qipongoing, admin, one_of(cov)) %>%
  na.omit()

mod_match <- matchit(qipongoing.binary ~ rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors,
                     method = "nearest",caliper = 0.2, data = temp_nomiss)

dta_m <- match.data(mod_match)
dta_m$qipongoing.binary <- as.factor(dta_m$qipongoing.binary)

summary(m1 <- glm.nb(rebel.violence.against.civilians.lead1 ~ qipongoing.binary + factor(admin), data=dta_m))
summary(m2 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.qipongoing + factor(admin), data=dta_m))
summary(m3 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.cost_qipongoing + factor(admin), data=dta_m))

stargazer(m1, m2, m3, 
          covariate.labels = c("QIP-binary", "QIPs-cost (logged)", "QIPs-number (logged)"), 
          add.lines = list(c("Fixed Effects", "District", "District", "District")),
          type="html", omit="admin", out="out.doc")


